Loading the packages and datasets

#Download packages
library(sf) 
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(sp)
library(tmap) 
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(rgeoda)
## Loading required package: digest
library(spdep)
## Loading required package: spData
## 
## Attaching package: 'spdep'
## The following object is masked from 'package:rgeoda':
## 
##     skater
#Set working directory
setwd("/Users/rebeccashi/Desktop/GGIS224/Final Project")

#load in chives
chives = st_read("chives-data-minusNU.geojson")
## Reading layer `chives-base' from data source 
##   `/Users/rebeccashi/Desktop/GGIS224/Final Project/chives-data-minusNU.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 801 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.94025 ymin: 41.64429 xmax: -87.52416 ymax: 42.02392
## Geodetic CRS:  WGS 84
#load in variables plant and human
plant = dplyr::select(chives,"trees_n")
heat = dplyr::select(chives,"heatisl")
human = dplyr::select(chives, "acs_population")
#summary
plant_summary = summary(plant$trees_n)
head(plant_summary)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##     0.000   874.000  1452.000  2135.905  2323.000 58668.000
heat_summary = summary(heat$heatisl)
head(heat_summary)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 23.00000 54.20000 71.30000 67.66205 82.90000 97.80000
human_summary = summary(human$acs_population)
head(human_summary)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##     0.000  2001.000  3077.000  3422.213  4500.000 19889.000
#plant_percentage
plant[is.na(plant)] <- 0
plant_sum = sum(plant$trees_n)
plant$plant_percentage = plant$trees_n/plant_sum

#human_percentage
human[is.na(human)] <- 0
human_sum = sum(human$acs_population)
human$human_percentage = human$acs_population/human_sum

#plant map
plant_map = tm_shape(chives) + tm_fill('trees_n', title = 'plant population',n=20, pal = "YlGn") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
          legend.title.size = 0.9,
          main.title = 'plant population, Chicago 2018', 
          main.title.size = 0.9)
print(plant_map)

#plant percentage map
plant_percentage_map = tm_shape(plant) + tm_fill('plant_percentage', title = 'percentage of plants',n=20,pal = "YlGn") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
          legend.title.size = 0.9,
          main.title = 'Percent population of plants, Chicago 2018', 
          main.title.size = 0.9)
print(plant_percentage_map)

#heat island index
heat_map = tm_shape(chives) + tm_fill('heatisl', title = 'heat island index',n=20, pal = "YlOrRd") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
          legend.title.size = 0.9,
          main.title = 'Heat island index, Chicago 2018', 
          main.title.size = 0.9)
print(heat_map)

#human map
human_map = tm_shape(chives) + tm_fill('acs_population', title = 'human population',n=20,pal = "RdPu") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
          legend.title.size = 0.9,
          main.title = 'human population, Chicago 2018', 
          main.title.size = 0.9)
print(human_map)

#human percentage map
human_percentage_map = tm_shape(human) + tm_fill('human_percentage', title = 'human percentage',n=20, pal = "RdPu") + tm_borders() + tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = 'right',
          legend.title.size = 0.9,
          main.title = 'Percentage of humans, Chicago 2018', 
          main.title.size = 0.9)
print(human_percentage_map)

#Spatial contiguity weight matrix PLANT
w.queen <- queen_weights(plant, order = 2, include_lower_order = TRUE)
lisa.queen <- local_moran(w.queen, plant["trees_n"], permutations = 999)
lisa_colors.queen <- lisa_colors(lisa.queen)
lisa_labels.queen <- lisa_labels(lisa.queen)
lisa_clusters.queen <- lisa_clusters(lisa.queen)
plot(st_geometry(plant), col=sapply(lisa_clusters.queen, function(x){return(lisa_colors.queen[[x+1]])}),  border = "#333333", lwd=0.2) 
title(main = "Plant population LISA (Queen)") 
legend('bottomleft', legend = lisa_labels.queen, fill = lisa_colors.queen, border = "#eeeeee")

lisa_p <- lisa_pvalues(lisa.queen)
p_labels <- c("Not significant", "p <= 0.05")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(plant), 
     col=sapply(lisa_p, function(x){
       if (x <= 0.05) return (p_colors[2])
       else return(p_colors[1])
       }), 
     border = "#333333", lwd=0.2)
title(main = "Plant Population LISA (Queen, P-Value)")
legend('bottomleft', legend = p_labels, fill = p_colors, border = "#eeeeee")

#Spatial contiguity weight matrix HEAT
w.queen <- queen_weights(heat, order = 2, include_lower_order = TRUE)
lisa.queen <- local_moran(w.queen, heat["heatisl"], permutations = 999)
lisa_colors.queen <- lisa_colors(lisa.queen)
lisa_labels.queen <- lisa_labels(lisa.queen)
lisa_clusters.queen <- lisa_clusters(lisa.queen)
plot(st_geometry(heat), col=sapply(lisa_clusters.queen, function(x){return(lisa_colors.queen[[x+1]])}),  border = "#333333", lwd=0.2) 
title(main = "heat island index LISA (Queen)") 
legend('bottomleft', legend = lisa_labels.queen, fill = lisa_colors.queen, border = "#eeeeee")

lisa_p <- lisa_pvalues(lisa.queen)
p_labels <- c("Not significant", "p <= 0.05")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(heat), 
     col=sapply(lisa_p, function(x){
       if (x <= 0.05) return (p_colors[2])
       else return(p_colors[1])
       }), 
     border = "#333333", lwd=0.2)
title(main = "heat island index LISA (Queen, P-Value)")
legend('bottomleft', legend = p_labels, fill = p_colors, border = "#eeeeee")

#Spatial contiguity weight matrix HEAT
w.queen <- queen_weights(human, order = 2, include_lower_order = TRUE)
lisa.queen <- local_moran(w.queen, human["acs_population"], permutations = 999)
lisa_colors.queen <- lisa_colors(lisa.queen)
lisa_labels.queen <- lisa_labels(lisa.queen)
lisa_clusters.queen <- lisa_clusters(lisa.queen)
plot(st_geometry(plant), col=sapply(lisa_clusters.queen, function(x){return(lisa_colors.queen[[x+1]])}),  border = "#333333", lwd=0.2) 
title(main = "pop density for humans LISA (Queen)") 
legend('bottomleft', legend = lisa_labels.queen, fill = lisa_colors.queen, border = "#eeeeee")

lisa_p <- lisa_pvalues(lisa.queen)
p_labels <- c("Not significant", "p <= 0.05")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(heat), 
     col=sapply(lisa_p, function(x){
       if (x <= 0.05) return (p_colors[2])
       else return(p_colors[1])
       }), 
     border = "#333333", lwd=0.2)
title(main = "pop density for humans  LISA (Queen, P-Value)")
legend('bottomleft', legend = p_labels, fill = p_colors, border = "#eeeeee")

#create dataset for plants and heat island index
hp = plant["plant_percentage"]
hp$heatisl = heat$heatisl
hp$pop = human$human_percentage

#correlation test with plant and heat island index
correlation = cor(hp$heatisl, hp$plant_percentage, method = "pearson")
head(correlation)
## [1] -0.08542622
#Overlay all the maps together
plant_bubble_map <- tm_shape(plant) + 
   tm_symbols(size = "plant_percentage", 
             col = "green", 
             alpha = 0.5)
heat_bubble_map <- tm_shape(heat) + 
   tm_symbols(size = "heatisl", 
             col = "red", 
             alpha = 0.5)
human_bubble_map <- tm_shape(human) + 
   tm_symbols(size = "human_percentage", 
             col = "blue", 
             alpha = 0.5)
all <- plant_bubble_map + heat_bubble_map + human_bubble_map
plantonly <- plant_bubble_map + heat_bubble_map
tmap_mode("plot")
## tmap mode set to plotting
print(plantonly)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

print(all)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

print(plant_bubble_map)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

print(heat_bubble_map)

print(human_bubble_map)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.